home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
FishMarket 1.0
/
FishMarket v1.0.iso
/
fishies
/
076-100
/
disk_085
/
plot6
/
plot6.c
< prev
next >
Wrap
C/C++ Source or Header
|
1992-05-06
|
6KB
|
239 lines
#include <stdio.h>
#include <math.h>
struct ent{
float ra,dc;
short mg,cl;
};
extern int MAX_X,MAX_Y;
/* 18.5cm/screen Y 26.2cm/screen X -> aspect ratio 1.4162 */
#define ASPECT 1.4162
/* MATH!
spherical to rectangular coordinates.
earth's axis is z axis
dc measured from positive (north) z axis
ra measured from x axis (we'll fix that later)
All stars are on a unit sphere.
x = sin(pi/2-dc)cos(ra)
y = sin(pi/2-dc)sin(ra)
z = cos(pi/2-dc)
rotation rule (about origin [ie about z axis])
a' = a cos(th) - b sin(th)
b' = a sin(th) + b cos(th)
rotate to see desired line of longitude(ra0).
ie line up y axis with chosen longitude (remember we were measuring from x)
rotate counter clockwise pi/2, rotate clockwise ra0 about z (ie pi/2-ra0)
x' = sin(pi/2-dc)cos(ra)cos(pi/2-ra0) - sin(pi/2-dc)sin(ra)sin(pi/2-ra0)
= sin(pi/2-dc)cos(ra+pi/2-ra0) [believe me!]
y' = sin(pi/2-dc)cos(ra)sin(pi/2-ra0) + sin(pi/2-dc)sin(ra)cos(pi/2-ra0)
= sin(pi/2-dc)sin(ra+pi/2-ra0)
z' = cos(pi/2-dc)
now rotate down from the north pole along your longitude to desired latitude
ie rotate down pi/2 then up by dc0
we are rotating around the x axis
x" = sin(pi/2-dc)cos(ra+pi/2-ra0)
y" = sin(pi/2-dc)sin(ra+pi/2-ra0)cos(pi/2-dc0)-cos(pi/2-dc)sin(pi/2-dc0)
= cos(dc)sin(ra+pi/2-ra0)cos(pi/2-dc0)-sin(dc)sin(pi/2-dc0)
z" = sin(pi/2-dc)sin(ra+pi/2-ra0)sin(pi/2-dc0)+cos(pi/2-dc)cos(pi/2-dc0)
= cos(dc)sin(ra+pi/2-ra0)sin(pi/2-dc0) + sin(dc)cos(pi/2-dc0)
I set dc0 to pi/2-dc0 and ra0 ot pi/2-ra0.
a = sin(dc) b = cos(dc)
c = sin(dc0) d = cos(dc0)
e = sin(ra + ra0) f = cos(ra + ra0)
x" = df;
y" = deb-ca;
z" = dea+cb;
Now! viewing transforms. I assume we are at (0,0,0) looking at ra0 x dc0.
We are looking down the z axis with y up.
our famous Ys = d/s Ye/Ze view box transform gives us:
(ie Y in screen coordinates from eye coordinates)
Xs = df/(dea+cb);
Ys = (deb-ca)/(dea+cb);
first we discard any stars behind us (ie z < 0)
if the star is more than half a screen to the side, discard it.
NB: we need not calculate f until the last moment.
divide x by the aspect ratio, then map to physical coordinates.
the rest is a subjective star plotting mechanism
*/
#define mabs(x) (((x)<0.0)?-(x):(x))
#define mfloor(x) (((x)<0.0)?-(int)-(x):(int)(x))
main(argc,argv)
int argc;
char **argv;
{
FILE *fp;
struct ent *ar,*star;
int t,i,num;
double a,b,c,d,e,f,de,ra0,de0,wi,xt,yt;
double HLF;
double z;
double TO_RADS;
if (argc <= 1) {
printf("usage plot <file-list>\n");
exit(0);
}
/* beware the ffp constant expression bug! */
TO_RADS = (TO_RADS=2.0) * 3.1415926535 / (TO_RADS=360.0);
HLF = 90.0 * TO_RADS;
printf("ra,de,wi [degrees]:");fflush(stdout);
scanf("%lf %lf %lf",&ra0,&de0,&wi);
ra0 = HLF-ra0*TO_RADS;
de0 = HLF-de0*TO_RADS;
wi = tan(mabs(wi)*TO_RADS);
window_init();
a = sin(de0); b = cos(de0);
while (argc-- > 1){
fp = fopen(argv[argc],"r");
if (fp == NULL){
printf("cant find %s\n",argv[argc]);
exit(0);
}
t = fread((char*)&num,4,1,fp);
if (t != 1) {
printf("didn't read file size!\n");
exit(0);
}
printf("reading %d records\n",num);
fflush(stdout);
ar = (struct ent*)malloc(num*sizeof(struct ent));
if (ar == NULL){
printf("couldn't malloc %d bytes\n",num*sizeof(struct ent));
exit(0);
}
t = fread((char*)ar,sizeof(struct ent),num,fp);
fclose(fp);
printf("read %d\n",t);
for (i = 0;i<t;i++) {
check_input();
star = ar+i;
c = sin(star->dc);
d = cos(star->dc);
e = sin(star->ra + ra0);
de = d * e;
z = wi*(de * a + c * b);
if (z <= 1.0e-16) continue; /* behind us, so skip it */
yt = (de * b - c * a) / z;
if (yt < -0.5 || yt > 0.5) continue;
f = cos(star->ra + ra0);
xt = d * f / ASPECT / z;
if (xt < -0.5 || xt > 0.5) continue;
plot_star((int)((xt + 0.5)*MAX_X),
(int)((yt + 0.5)*MAX_Y),star->mg,star->cl);
}
printf("done\n");
fflush(stdout);
free (ar);
}
while (1){
wt();
check_input();
}
}
plot_star(x,y,mg,cl)
int x,y;
short mg,cl;
{
int c,colour,m;
m = (int)mfloor(((double) mg)/100.0+0.49);
c = (int)mfloor(cl/60.0+0.49) + 1;
if (c < 0) c = 0;
else if (c > 3) c = 3;
c = 1 + 3 * c + 2;/*dimmest of range*/
colour = c-1;/* medium */
switch(m){
case 0: point(x-3,y,colour);
point(x-4,y,colour);
point(x+3,y,colour);
point(x+4,y,colour);
point(x+5,y,colour);
point(x-5,y,colour);
colour = c-2;
point(x+2,y+1,colour);
point(x+2,y-1,colour);
point(x-2,y+1,colour);
point(x-2,y-1,colour);
case 1: point(x,y+2,colour);
point(x,y-2,colour);
colour = c-2;
case 2: /* in bright colour */
point(x+3,y,colour);
point(x-3,y,colour);
point(x+2,y,colour);
point(x-2,y,colour);
colour = c-2;
point(x+1,y+1,colour);
point(x+1,y-1,colour);
point(x-1,y+1,colour);
point(x-1,y-1,colour);
/* in brightest colour */
case 3: /* in dimmest colour */
colour = c-2;
/* brighten colour */
case 4: /* in dimmest colour*/
point(x,y+1,colour);
point(x,y-1,colour);
point(x+1,y,colour);
point(x-1,y,colour);
point(x+2,y,colour);
point(x-2,y,colour);
case 5: colour = c-2;/* set to brightest */
break;
case 6: colour = c-1;
break;
default: colour = c;/* in dimmest colour */
}
point(x,y,colour);
}